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ABSTRACT 

We present the results of both global cylindrical disc simulations and local shearing 
box simulations of protoplanets interacting with a disc undergoing MHD turbulence 
with zero net flux magnetic fields. We investigate the nature of the disc response and 
conditions for gap formation. This issue is an important one for determining the type 
and nature of the migration of the protoplanet, with the presence of a deep gap being 
believed to enable slower migration. 

For both types of simulation we find a common pattern of behaviour for which 
the main parameter determining the nature of the response is M P R 3 / (M„H 3 ) 1 with 
M P ,M*,R, and H being the protoplanet mass, the central mass, the orbital radius 
and the disc semi-thickness respectively. We find that as M p R 3 / (A1*H 3 ) is increased 
to ~ 0.1 the presence of the protoplanet is first indicated by the appearance of the 
well known trailing wake which, although it may appear to be erratic on account of 
the turbulence, appears to be well defined. Once M p R 3 /(M*H 3 ) exceeds a number 
around unity a gap starts to develop inside which the magnetic energy density tends 
to be concentrated in the high density wakes. This condition for gap formation can 
be understood from simple dimensional considerations of the conditions for nonlinear- 
ity and the balance of angular momentum transport due to Maxwell and Reynolds' 
stresses with that due to tidal torques applied to the parameters of our simulations. 

An important result is that the basic flow morphology in the vicinity of the pro- 
toplanet is very similar in both the local and global simulations. This indicates that, 
regardless of potentially unwanted effects arising from the periodic boundary con- 
ditions, local shearing box simulations, which are computationally less demanding, 
capture much of the physics of disc-planet interactions. Thus they may provide a use- 
ful tool for studying the local interaction between forming protoplanets and turbulent, 
protostellar discs. 

Key words: accretion, accretion disks — MHD, instabilities, turbulence - planetary 
systems: formation, protoplanctary discs 



1 INTRODUCTION 

The recent and ongoing discovery of extrasolar giant plan- 
ets has stimulated much current investigation of theories 
of planet formation (e.g. Mayor & Queloz 1995; Marcy, 
Cochran, & Mayor 1999; Vogt et al. 2002). At the present 
time giant planets are believed to originate in protostellar 
discs. Possible formation scenarios are either through direct 
gravitational fragmentation of a young protostellar disc (e.g. 
Boss 2001) or through the accumulation of a solid core which 
then undergoes mass build up through gas accretion once the 
core mass reaches ~ 15 Earth masses (e.g. Bodenheimer & 
Pollack 1986; Pollack et al. 1996). In either scenario giant 



planets are likely to form at significantly larger radii than 
those observed, requiring a migration mechanism to bring 
them closer to the central star. Disc-protoplanet interaction 
provides a natural mechanism. 

In the standard picture a protoplanet exerts torques on 
a protostellar disc through the excitation of spiral density 
waves (e.g. Goldreich & Tremaine 1979). These waves carry 
an associated angular momentum flux which is deposited in 
the disc where the waves are damped. This process results in 
a negative torque acting on the protoplanet from the outer 
disc and a positive torque acting on it from the disc interior 
to its orbit, with the outer disc torque normally dominating. 
For low mass planets that induce linear perturbations in the 
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disc, the resulting migration is usually referred to as type I 
migration. 

Previous work on the interaction between a protoplanet 
and a laminar disc with Navier Stokes viscosity (Papaloizou 
& Lin 1984; Lin & Papaloizou 1993) indicates that a suf- 
ficiently massive protoplanet can open up an annular gap 
in the disc centred on its orbital radius. For typical pro- 
tostellar disc models the protoplanet needs to be approxi- 
mately a Jovian mass for gap formation to occur. Recent 
simulations (Bryden et al. 1999; Kley 1999; Lubow, Seib- 
ert, & Artymowicz 1999; D'Angelo, Henning, & Kley 2002) 
examined the formation of gaps by giant protoplanets, and 
also estimated the maximal gas accretion rate onto them. 
The orbital evolution of a Jovian mass protoplanet embed- 
ded in a standard laminar viscous protostellar disc model 
was studied by Nelson et al. (2000). They found that that 
gap formation and accretion of the inner disc by the cen- 
tral mass led to the formation of a low density inner cavity 
in which the planet orbits. Interaction with the outer disc 
resulted in inward migration on a time scale of a few x 10 5 
yr. Here the viscous evolution of the disc is responsible for 
driving the inward migration of the giant planet. This mode 
of migration is normally referred to as type II migration. 
Recent work examining gap formation, mass accretion, and 
migration torques in three dimensions has been presented by 
D'Angelo, Kley, & Henning (2003) and Bate et al. (2003). 
The disc models in these studies all adopted an anomalous 
disc viscosity modelled through the Navier-Stokes equations 
without consideration of its origin. 

The most likely origin of the viscosity is through MHD 
turbulence resulting from the magnetorotational instability 
(MRI) (Balbus & Hawley 1991) and it has recently become 
possible through improvements in computational resources 
to simulate discs in which this underlying mechanism re- 
sponsible for angular momentum transport is explicitly cal- 
culated. This is necessary because the turbulent fluctuations 
do not necessarily result in transport phenomena that can 
be modelled with the Navier Stokes equation. 

To this end Papaloizou & Nelson (2003) and Nelson & 
Papaloizou (2003a) (hereafter papers I and II) developed 
models of turbulent protostellar accretion discs and consid- 
ered the interaction with a giant protoplanet of 5 Jupiter 
masses. The large mass was chosen to increase the scale of 
the interaction so reducing the computational resources re- 
quired. This protoplanet was massive enough to maintain 
a deep gap separating the inner and outer disc and exert 
torques characteristic of type II migration. A study of gap 
formation in turbulent discs has also been presented by Win- 
ters, Balbus, & Hawley (2003a) 

It is the purpose of this paper to extend our previous 
work to a wider range of protoplanetary masses and disc 
parameters. In particular we study the nature of the disc 
protoplanet interaction as the protoplanet mass increases 
and investigate the condition for gap formation. 

To do this we utilise both global and local simulations. 
The global simulations are more realistic but too computa- 
tionally demanding for significant extended parameter sur- 
veys. On the other hand local shearing box simulations being 
less computationally demanding enable a wider parameter 
survey at higher resolution. Of course because of the high 
degree of symmetry they do not allow migration rates to be 



estimated but most other features of the global simulations 
are reproduced. 

In this paper we focus on the morphology of the disc 
response as a function of the protoplanet mass leaving the 
analysis of the disc-protoplanet torques resulting from the 
interaction to a companion paper (Nelson & Papaloizou 
2003b - hereafter paper IV). 
The plan of this paper is as follows: 

In §2 we describe both the initial global disc models and the 
shearing box models used for the simulations. 
In §3 we discuss the nature of the waves that can propa- 
gate and how their excitation leads to the existence of the 
prominent wake that is frequently seen in the simulations of 
protoplanet disc interaction even when the disc is turbulent. 
In §4 we describe the numerical procedure used. 
In §5 we present the results of the local and global simula- 
tions. Finally, we summarize our results in §6 . 



2 SET UP OF INITIAL MODELS 

In this paper we have performed both global and local simu- 
lations of protoplanets interacting with model discs in which 
MHD turbulence driven by the MRI operates. We consider 
simulations for which the magnetic field has zero net flux and 
therefore an internally generated dynamo is maintained. We 
determine the disc response as a function of the protoplanet 
mass and as far as is possible study its characteristics for the 
different kinds of model focusing on the conditions that de- 
termine whether the response is weak (linear) or strong (non 
linear). We now describe the model set up for the global and 
local simulations. 



2.1 Global Simulations 

The governing equations for MHD written in a frame ro- 
tating with arbitrary uniform angular velocity f2 p k with k 
being the unit vector along the rotation axis assumed to be 
in the vertical (z) direction are the continuity equation: 



the equation of motion 
<9v 



dt 



+ v ■ Vv + 20,kxv = - 



(1) 



Vp „^ (V x B) x B , . 
— - V$ + ^ —!■ , 2 

p Anp 



and the induction equation 



dB 

dt 



= V x (v x B). 



(3) 



where v, P, p, B and $ denote the fluid velocity, pressure, 
density and magnetic field respectively. The potential $ = 
®rot + contains contributions due to gravity, 3>g and the 
centrifugal potential $ ro t = — (l/2)fl P |k x r| 2 . 

The gravitational contribution is taken to be due to a 
central mass M» and planet with mass M p . Thus in cylindri- 
cal coordinates (r, 0, z) with the planet located at (r p , cj> p , 0) 
and the star located at (r», 0», 0) $g = $c + $ P , where 



$V = - 



GM, 



and <!>„ 



GM V 



y/r 2 + r 2 — 2rr p cos(0 — <j> p ) + b 2 
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Here, as in papers I and II, for reasons of computational ex- 
pediency, we have neglected the dependence of the gravita- 
tional potentials on z along with the vertical stratification of 
the disc. Thus the global simulations are of cylindrical discs 
(e.g. Armitage 1998, 2001; Hawley 2000, 2001; Steinacker & 
Papaloizou 2002). To model the effects of the reduction of 
the planet potential with vertical height, we have incorpo- 
rated a softening length b in its potential. 
We use a locally isothermal equation of state in the form 

P = p-c(r) 2 , (5) 

with c denoting the sound speed which is specified as a fixed 
function of the cylindrical radius r. The above gives the basic 
equations used for the global simulations. 

2.1.1 Initial and Boundary Conditions 

The global simulations presented here all use the same un- 
derlying turbulent disc model. This model has a constant 
aspect ratio H/r = c(r)/(rQ) = 0.07, where Q is the disc 
angular velocity as measured in the inertial frame. The inner 
radial boundary of the computational domain is at Ri n — 1 
and the outer boundary is at R ou t = 8. The boundary condi- 
tions employed are very similar to those described in paper 
I. Regions of the disc in the close vicinity of the inner and 
outer boundaries were given non-Keplerian angular velocity 
profiles that are stable to the MRI, and which have large 
values of the density in order to maintain radial hydrostatic 
equilibrium. These regions act as buffer zones that prevent 
the penetration of magnetic field to the radial boundaries, 
thus maintaining zero-net flux in the computational domain. 
The inner buffer zone runs from r = 1.14 to r = Ri n , and 
the outer buffer zone runs from r — 7.0 to r = R ou t- 

The disc was initiated with a zero net flux toroidal mag- 
netic field in a finite annulus in the disc between the radii 
2 < r < 6. The initial field varied sinusoidally over four com- 
plete cycles within this radial domain, and the initial ratio 
of the volume integrated magnetic energy to the volume in- 
tegrated pressure was (1//3) = 0.032. The initial disc model 
had azimuthal domain running between 4> = to 7r/4, and 
vertical domain from z min = —0.21 to z max — 0.21. Peri- 
odic boundary conditions were employed in the vertical and 
azimuthal directions. 

The unperturbed disc model was evolved until the MRI 
produced a turbulent disc model with a statistically steady 
turbulence. This model had a volume averaged value of the 
ratio of magnetic energy to mean pressure of ~ 0.011 and 
a volume averaged value of the Shakura-Sunyaev stress pa- 
rameter a ~ 7 x 10~ 3 . In order to generate the full 2n disc 
models used to study the interaction between turbulent discs 
and protoplanets, eight identical copies of this relaxed disc 
sector were joined together. The model described in table 2 
with azimuthal domain ty/2 (run G5) contained only two 
copies of the 7r/4 sector, and periodic boundary conditions 
were employed in the azimuthal domain. 

Gravitational softening of the protoplanet was em- 
ployed in all simulations. In runs G1-G4 described in ta- 
ble 2, the softening parameter b — 0.331/. For run G5, this 
was reduced to b = 0.1H. For runs Gl, G2, and G3, the 
planet was placed at a radius such that \r p — r*| =3 and 
at azimuthal location cj> p = n in the disc. For runs G4 it 
was placed such that (\r p — r*\,(f> p ) = (2.5, tv) and for run 



G5 it was placed at (\r p — r*\,4> p ) = (2.5, 7r/4). All simu- 
lations were performed in a rotating frame whose angular 
velocity equalled that of the circular Keplerian planetary 
orbit. The unit of time adopted when discussing the results 
is Q^ 1 evaluated at r — 1, the inner boundary of the com- 
putational domain. We note that the orbital period at this 
radius is P(r = 1) = 2ir, the orbital period of a planet at 
r p = 2.5 is P(r — 2.5) ~ 24, and the orbital period of a 
planet at r p = 3 is P(r = 3) ~ 32. 

In order to explore the disc-planet interaction over a 
wider range of parameter space and with higher resolution 
in the vicinity of the planet, we have adopted a local shearing 
box approximation. This we describe below. 

2.2 Local Shearing Box Simulations 

In the shearing box limit (Goldreich & Lynden-Bell 1965) 
we consider a small Cartesian box centred on the point at 
which the Keplerian angular velocity is £l p . We define local 
Cartesian coordinates (a;, y, z) with associated unit vectors 
(i, j, k). The direction i is taken to be outward along the line 
joining the origin to the central object while k points in the 
vertical direction. The direction j is along the unperturbed 
rotational shear flow. 

The equation of motion may then be written as 

^ + vW+2fi p kxv-3fi^i = -gP-yg p+ ( V x B ) xB ( 6 ) 
at p 4:Tvp 

The induction and continuity equations retain the same 
form as equations (3) and (1) but are of course expressed in 
component form in the local box centred Cartesian coordi- 
nates. We recall that the term oc x in equation (6) is derived 
from a first order Taylor expansion about x — of the com- 
bination of the gravitational acceleration due to the central 
mass and the centrifugal acceleration (Goldreich & Lynden- 
Bell 1965). However, as for the global disc simulations, wc 
have neglected the z dependence of the gravitational poten- 
tial and therefore vertical stratification. Thus in this respect 
the simulations are of unstratified boxes of the type consid- 
ered by Hawley, Gammie & Balbus (1995). 

Here, when introduced, the planet is located in the cen- 
tre of the box, thus 



In the steady state box with no planet or magnetic field, 
the equilibrium velocity is due to Keplerian shear and thus 

v = (0,-3fi p a:/2,0). (8) 

We adopt an isothermal equation of state P = pc 2 , 
with c being the sound speed which is taken to be constant 
throughout the box. 

2.2.1 Initial and Boundary Conditions 

The shearing box is presumed to represent a local patch of a 
differentially rotating disc. Thus the appropriate boundary 
conditions on the bounding faces y = constant = ±Y and 
z — constant — ±Z derive from the requirement of period- 
icity in the local Cartesian coordinate directions normal to 
the boundaries. On the boundary faces x = constant = ±X, 
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the boundary requirement is for periodicity in local shear- 
ing coordinates. Thus for any state variable F(x,y,z), the 
condition is that F(X,y,z) = F(~X,y - 3Q p Xt,z). This 
means that information on one radial boundary face is com- 
municated to the other boundary face at a location in the 
azimuthal coordinate y shifted by the distance the faces have 
sheared apart since the start of the simulation, t, (see Haw- 
ley, Gammie & Balbus 1995). 

For all models adopting a shearing box, the softening 
parameter was taken as b — 0.3H. The potential was flat- 
tened (made to attain a constant value in a continuous man- 
ner) at large distances from the protoplanet in order to sat- 
isfy the periodic boundary conditions. In simulations BaO- 
Ba4 (described in table 1) this was done outside the cylinder 
of radius 3H centred on the planet. In simulations Bbl-Bb4 
this was done only close to the boundary. 

The simulations Ba0-Ba4 occupied a total width 8-ff in 
x and 4ttH in y. These boxes were found to be large enough 
to accomodate the scale of the disc response to forcing by the 
protoplanet and enable gap formation. But note that inter- 
ference from protoplanets in neighbouring boxes implied by 
the periodic boundary conditions is not eliminated entirely. 
In the small protoplanet mass regime waves propagate be- 
tween boxes and in the nonlinear regime torques exerted by 
protoplanets in neighbouring boxes may influence the gap 
formation process. However, comparison between local and 
global runs indicates that such effects are not major and 
do not prevent the box runs from capturing the essential 
physics. 

To examine the effect of varying box size we ran simu- 
lations Bbl-Bb4 which were of the same resolution but had 
twice the extent in y. 

All of these models were initiated by imposing a vertical 
field, with zero net flux through the box, varying sinusoidally 
in x with a wavelength of H. The amplitude was such that 
the initial value of the ratio of the total magnetic energy to 
volume integrated pressure was 0.0025. The initial velocity 
in the x direction at each grid point was chosen to be the 
product of a random number between —1.0 and 1.0 and 0.1c. 
The initial velocity in the y-direction is given by equation 8, 
and the z component of velocity was initialized to zero. The 
unit of time for all simulations was taken to be Q.p 1 . In these 
units the orbital period at the centre of the box is 2tv. 

2.3 Dimensionless Variables and Scaling with 
Planet Mass 

It is helpful to write the governing equations for the shear- 
ing box in dimensionless form. In this way the dependence 
of the outcomes of the simulations on planet mass and disc 
parameters is made clearer (see also Korycansky & Pa- 
paloizou 1996). To do this we adopt dimensionless coordi- 
nates (x',y',z') — (x/ H,y /H, z/H), where the length scale 
H — c/flp would correspond to the disc scale height were 
it vertically stratified and the dimensionless time t' = Q. v t. 
We also note that the Keplerian relation Q 2 = GM*/R 3 , 
may be used to relate the rotation rate of the centre of the 
box to a putative central mass M, and orbital radius R. 
We now introduce the dimensionless velocity, sound speed 
and planet mass through v' = v/(fl p H), d — c/(fl p H), and 
M' p = M P R 3 /(M,H 3 ) respectively. As we do not include the 
self-gravity of the disc in calculating its response, the mag- 



nitude of the disc density may be arbitrarily scaled. The 
models simulated here all start with a uniform density po 
which, together with H, may be used to specify the mass 
scaling. Then we define a dimensionless density and mag- 
netic field through p' — pop and B' = R/(£l p Hy/4irp ). The 
equation of motion may be written in the above dimension- 
less variables in the form 

4— + V-X7V + 2kxv' - 3x'i = 
of 

+ (9) 

where 

/ M' , s 

K = ; (10) 

^Jx' 2 + y' 2 + b' 2 

with dimensionless softening parameter b' — b/H and of 
course the spatial derivatives are with respect to the dimen- 
sionless coordinates. Similarly the dimensionless induction 
equation becomes 

EEL = V x (v 7 x B'). (11) 
and the dimensionless continuity equation is 
g^+V(.p'v')=0. (12) 

We remark that in dimensionless coordinates the 
boundaries of the box are x' = ±X/H, y' — ±Y/H, 
z' — ±Z/H. Also if the magnetic field has zero net flux 
and is dynamo generated and thus a spontaneous prod- 
uct of the simulation, given that the only parameter oc- 
curring in equations (9 - 12 ) is the dimensionless mass 
M' v = M P R 3 /(M*H 3 ), we should be able to consider the 
dependence of the time averaged outcome of a simulation to 
be only on X/H, Y/H, Z/H, and M' p = M P R 3 /(M*H 3 ). 

For a planet embedded in a box with natural verti- 
cal scale H with sufficiently large X/H to incorporate the 
tidal interaction without interference from other boxes cen- 
tred on different radii, the only parameters determining 
long term averaged properties of a simulation should be 
Y/H and M' v . Thus for such a quantity Q, we should have 
Q = f(Y/H,M p R 3 /(M t H 3 )). 

Note that the above discussion implies that for a box 
of fixed dimension, the only distinguishing parameter is 
M' v — MpR 3 /(M f H 3 ). This parameter may also be inter- 
preted as the cube of the ratio of the Hill radius to disc 
scale height and the condition that M' v > 1 leads to the 
thermal condition of Lin & Papaloizou (1993) that the Hill 
radius should exceed the disc scale height for gap formation 
to occur. From Korycansky & Papaloizou (1996) this con- 
dition is also required in order that the perturbation due to 
the protoplanet be non linear. 

The above discussion indicates that were we to vary the 
box size in the azimuthal (y) direction, the condition for gap 
formation should be of the form M' v > f g (Y/H). As gap 
formation by the same planet in a longer box is expected 
to be more difficult, we expect f g to be a non decreasing 
function of its argument. In order to correspond to a full 
circle of radius R, we require that Y = -kR. 

It is perhaps suggestive to point out that if f g {irR/H) = 
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Model z domain a; domain y domain ti (M P R 3 /(M t H 3 ) n z n x ra 



BaO 


(- 


-if/2, if/2) 


{-AH, AH) 


(-2-kH, 2txH) 


Bal 


(- 


-H/2, if/2) 


(—AH, AH) 


(-2itH, 2-kH) 


Ba2 


(- 


-H/2, H/2) 


(—AH, AH) 


(-2-kH, 2-kH) 


Ba3 


(" 


-H/2, H/2) 


(—AH, AH) 


(-2itH, 2-kH) 


Ba4 


(" 


-H/2, H/2) 


(—AH, AH) 


(-2nH, 2-kH) 


Bbl 


(" 


-H/2, H/2) 


(—AH, AH) 


(-AttH, AwH) 


Bb2 


(- 


-H/2, H/2) 


(—AH, AH) 


(-4-kH, AttH) 


Bb3 


(- 


-H/2, H/2) 


(—AH, AH) 


(-AttH, AttH) 


Bb4 


(" 


-H/2, H/2) 


(—AH, AH) 


(-AttH, AttH) 



0.0 


637 


0.0 


35 


261 


200 


354 


621 


0.1 


35 


261 


200 


354 


644 


0.3 


35 


261 


200 


354 


619 


1.0 


35 


261 


200 


354 


650 


2.0 


35 


261 


200 


216 


392 


0.2 


35 


261 


396 


216 


431 


0.5 


35 


261 


396 


216 


411 


1.0 


35 


261 


396 


216 


485 


2.0 


35 


261 


396 



Table 1. Parameters of the shearing box simulations: The first column gives the simulation label, the second, third and fourth give the 
extent of the coordinate domains considered. The first nine simulations were performed in shearing boxes with the x, y, and z domains 
referring to the Cartesian domains. The last two simulations were global so that in these cases the x, y, and z domains refer to the r, <j>, and 
2 domains for the cylindrical coordinate system used. The fifth and sixth column give the start and end times if the simulation measured 
in dimcnsionlcss units. Simulations with protoplanets were started by inserting the protoplanet into the turbulent model generated from 
the simulation BaO at time indicated. The seventh column gives the value of (G M p / (Q, 2 H 3 ) . The eighth, ninth and tenth columns give 
the number of computational grid points in the x,y, and z coordinates respectively. These numbers include any ghost zones used to 
handle boundary conditions. 



40a(R/H), the condition for gap formation indicated above 
becomes M p /M* > 40aH 2 /(R 2 ) which is of the same form 
as the viscous criterion of Lin & Papaloizou (1993) which 
was obtained by requiring angular momentum transport due 
to excited waves to exceed that due to viscosity. However, 
it is important to note that if we introduce such an a pa- 
rameter here, it cannot correspond to an imposed viscosity 
as in Lin & Papaloizou (1993). Instead it must be related to 
the turbulent transport arising from MHD turbulence. We 
take it to correspond to a time average of the volume av- 
eraged stress {B x B y /(4n)) expressed in units of the volume 
averaged pressure (P) , thus corresponding to the Shakura & 
Sunyaev (1973) parameterization. This procedure is reason- 
ably well defined in the parts of the disc not strongly per- 
turbed by a protoplanet (eg. Paper II) . For simulations with 
zero net magnetic flux a~5x 1CP 3 but it takes on larger 
values when a net poloidal or toroidal flux is imposed (eg. 
Hawley, Gammie & Balbus 1996; Steinacker & Papaloizou 
2002). 

Adopting a — 5 x 10~ 3 , suggests f g (Y/H) ~ 
0.07(Y/H). We recall that should 0.07(Y/H) be less than 
unity, the requirement of nonlinear perturbation by the pro- 
toplanet in order for gap formation to occur is M p /M* > 
C t (H/R) 3 , where Ct is a constant of order unity (Korycan- 
sky & Papaloizou 1996). Both these conditions can be com- 
bined to give the requirement that 

M P R S /(M^H 3 ) > max(C t , 0.07(Y/H)) (13) 

In this way we may see how both the thermal and viscous 
conditions of Lin & Papaloizou (1993) may be related to the 
simulations. 



3 WAKES AND CHARACTERISTICS 

A well known phenomenon found in simulations of proto- 
planets interacting with laminar discs is a prominent wake 
that is observed to trail away from the protoplanet eventu- 



ally becoming a spiral wave propagating away from it (e.g. 
Kley 1999; Nelson et al 2000; Lubow, Seibert & Artymowicz 
1999; Ogilvie & Lubow 2002). Such structures, though dis- 
torted, are also seen in our simulations of discs with MHD 
turbulence in which there are strong and highly time vari- 
able fluctuations, even for weakly perturbing protoplanets. 
In view of their importance and persistence we here give a 
brief discussion of their interpretation as being defined by a 
characteristic ray emanating from the protoplanet that sep- 
arates the flow shearing past it into parts causally connected 
and not causally connected by waves. The wakes thus natu- 
rally correspond to the location a strong wave sent out from 
the protoplanet into the flow shearing past it. 



3.1 The Characteristic Rays 

The set of basic equations defined by (1), (2) and (3) are 
well known to form a hyperbolic set which can be written 
in the standard form (Courant & Hilbert 1962) 

E^gW^O (* = l,2,.... > »). (14) 

Here the set of equations (1), (2) and (3) are arranged into 
a set of n simultaneous first order PDEs for the n compo- 
nents Ui of U which are in turn composed of p together 
with the components of B and v. The Xj (j = 1 — 4) are the 
three spatial coordinates and the time t respectively. The 
characteristic rays, applicable to both the global and local 
formulation may be found by first obtaining the local dis- 
persion relation associated with (14). To do this one sets 
Ui = Uio exp(ikjXj — iiot) under the assumption that the 
magnitudes of the components of the wave vector ki(l — 
1, 2, 3) = k = (k x ,ky,k z ) and the wave frequency |a>| arc 
very large. The dispersion relation is then given by the van- 
ishing of a determinant of coefficients \A^kj — A^u>\ — 0. 
As is very well known (eg. Campbell 1997) this dispersion 
relation leads to three distinct wave modes. The first are the 
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Model <j> domain H/r M p /M* (M P R 3 /(M„H 3 ) n r n z 



Gl 


2tt 


0.07 


1 x 10~ 5 


0.03 


450 


1092 


40 


G2 


2tt 


0.07 


3 x 10~ 5 


0.09 


450 


1092 


40 


G3 


2tt 


0.07 


1 x 10~ 4 


0.30 


450 


1092 


40 


G4 


2tt 


0.07 


3 x 10~ 3 


8.75 


450 


1092 


40 


G5 


tt/2 


0.07 


3 x 10~ 3 


8.75 


450 


276 


40 



Table 2. Parameters of the global simulations: The first column gives the simulation label, the second gives the extent of the azimuthal 
domain, the third gives the H/r value of the disc, and the fourth gives the protoplanet-star mas ratio. The fifth gives the ratio of M p /M* 
to (H/r) 3 . The sixth, seventh, and eighth columns describe the number of grid cells used in each coordinate direction. 



Alfven waves propagating in opposite directions which have 
u> = k • v ± k • va, with Vi = «aB/|B|, where va is the 
Alfven speed being such that v\ = |B| 2 /(47rp). The other 
modes are fast and slow magnetosonic waves which satisfy 



k\c 2 +v A )Jk^ + v A ) 



2fc 2 c 2 (k- v A ) 2 



(u;+k-v) 2 = 

Here k = |k| and we note that the frequency, u> occurs added 
to the quantity k • v in these relations. This is due to the 
Doppler shift occurring on account of the motion of the fluid. 

For each of these waves we may specify u> = w(k, r, t). 
Then the characteristic rays follow from integrating the 
equations of geometric optics (e.g. Whitham 1999) 

(16) 



.(15) 



dr du) dk. 
~r = and —— 

dt 9k dt 



dr 



To be specific, we here concentrate on the fast mag- 
netosonic mode which in our case leads to the most rapid 
communication. Further, as for the most part the magnetic 
energy is on the order of one percent of the thermal energy, 
we shall neglect the field so that the mode becomes an or- 
dinary sound wave such that 



-k • v + ck. 



(17) 



We consider the general situation to be one for which 
the flow streams past the protoplanet from y > to y < 
for x > and with the flow moving in the opposite sense 
for x < 0. Considering without loss of generality the case 
y > 0, when the flow streams past supersonically, if we ne- 
glect dependence on z, it is divided into two parts by the 
characteristic originating close to the planet from the point 
with y = and the smallest possible positive value of x at 
which the flow speed is equal to the wave speed. If the planet 
had no gravity, the fluid that is upstream of this character- 
istic would be unaware of the presence of the planet. 

The fluid that first knows about the planet, namely that 
on the characteristic is a natural channel for waves to prop- 
agate that are excited by the planet which is why this region 
is so prominent in the simulations. Note that the concept of 
characteristic rays is robust in that it can apply to a general 
time dependent fluid with turbulence. Also with z depen- 
dence restored we would expect the fluid to be divided by 
a surface made up from many rays rather than a single ray 
which is of the same form independent of z. This would lead 
to some blurring of the wake. 

A characteristic ray is defined by its point of origin (here 
specified above) and the value of k. However, note that only 



the ratio of the components and not the magnitude matters. 
The natural choice is k z = 0, and k y /k x such that to — 
0. This is because we expect a zero frequency wave to be 
excited by the planet in a frame in which that is at rest. 
However, note that in a turbulent time dependent situation 
lo would not be constant. 

The simplest situation is one when the flow is time in- 
dependent and the velocity is given by (8). In that case it 
follows from equations (16) that both u> = and k y are 
constant. 

It further follows from equations (16) that the charac- 
teristic rays then satisfy 



dx 



= -y/v* /<?-!. 



(18) 



Using equation (8) for v y and integrating we obtain for 
the characteristic ray for x > xo 

V ^ X ° = _ 2a " v/ ^ 5 ^° ~ 1+ \ ln (v^V^o ~ l + x/xoj ,(19) 

with xo = 2H/3. To complete the picture we remark that 
the characteristic for x < 0, can be found from the condition 
y(-x) = -y(x). 

We comment that the computational domains in all of 
our simulations are extensive enough to contain the char- 
acteristic associated with the waves emitted by an exciting 
object, be it turbulent fluctuations or a protoplanet (i.e. x 
exceeds 2/3H by a significant margin). 



4 NUMERICAL PROCEDURE 

The numerical scheme that we employ is based on a spa- 
tially second-order accurate method that computes the ad- 
vection using the monotonic transport algorithm (Van Leer 
1977). The MHD section of the code uses the method of 
characteristics constrained transport (MOCCT) as outlined 
in Hawley & Stone (1995) and implemented in the ZEUS 
code. The code has been developed from a version of NIR- 
VANA originally written by U. Ziegler (see SP, Ziegler & 
Riidiger (2000), and references therein) 
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Figure 1. The ratio of the total magnetic energy to the volume 
integrated pressure, (1//3), is plotted as a function of time for 
models BaO, Bal, Ba2, Ba3, and Ba4. The plots begin to separate 
after the protoplanet is inserted at time 354.. At time 500 the 
curves correspond to simulations Ba4, Bal, Ba2,Ba3 and BaO 
with decreasing values of (1//3) respectively. Although there are 
large fluctuations, the maximum variation at any time is factor 
of 3. This is significantly smaller than the density depression in 
the deepest gap (see figure 10). 




Figure 2. The volume averaged stress parameter (a) is plotted as 
a function of time for model BaO which has no protoplanet. The 
noisy full curve contains contributions from both the magnetic 
and Reynolds' stresses while the dotted curve gives the magnetic 
contribution. The dashed curve gives the running time average 
of the total stress. The large initial values are due to the initial 
strong instability which produces a channel phase. However, af- 
ter about four orbits at the centre of the box, the running time 
average is within a factor of two of the final value. 



Figure 3. Density contour plot in a (x, y) plane for simulation 
BaO which has no embedded protoplanet. The extensions are 8H 
in the x direction and 4irH in the y direction. The plot applies 
near the end of the simulation. The range of variation in the den- 
sity is about a factor of three. This is typical of what is obtained 
in any (x,y) plane once a statistically steady state turbulence is 
set up. Note the trailing elongated structures. 



4.1 Vertically and Horizontally Averaged Stresses 
and Angular Momentum Transport 

In order to describe average properties of the turbulent mod- 
els, following papers I and II we use quantities that are ver- 
tically and azimuthally averaged over the (y, z) domain for 
the shearing boxes or the ((f>, z) domain for the global sim- 
ulations (e.g. Hawley 2000). Thus we may consider y = <f> 
and x = r in these cases. Sometimes an additional time av- 
erage may be adopted. The vertical and azimuthal average 
of some quantity Q is then defined through 



7T / 



zdy 



J pdzdy 

The disc surface density is given by 



pdzdy. 



(20) 



(21) 



The vertically and azimuthally averaged Maxwell and 
Reynolds stresses, are given by: 



Tm = 2YE ' BxBy 



4np 
and 

Tr £ — 2V£8v x Sv y . 



(22) 



(23) 



Here the velocity fluctuations 5v x and 8v y are defined 
through, 



Sv x 



SVy = Vy~Vy 



(24) 
(25) 



The horizontally and vertically averaged Shakura & Sunyaev 
(1973) a stress parameter appropriate to the total stress is 
given by 
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Figure 5. The volume averaged total stress parameter (a) is 
plotted as a function of time for models Bal, Ba2, Ba3, and Ba4. 
The plots begin from the moment the protoplanet is inserted. The 
noisy upper curves contain contributions from both the magnetic 
and Reynolds' stresses. They can be identified at early times as for 
Ba4 (uppermost full curve), Ba3 (uppermost dotted curve), Ba2 
(next highest full curve), Bal (next highest dotted curve). The 
lower and much less noisy curves give the magnetic contribution. 
These are similar until they separate at about time 500. At this 
time they may be identified with increasing values of the stress 
parameter as associated with Bal,Ba3, Ba2 and Ba4 respectively. 
The values of the total (a) show a sharp increase when the proto- 
planet is inserted in the cases with higher mass protoplanets on 
account of the excitation of density waves which provide angular 
momentum transport. 




5 NUMERICAL RESULTS 



Figure 4. Contours of magnetic energy density corresponding 
to the plot in figure 3. The upper plot corresponds to the (x,y) 
plane illustrated in figure 3 while the lower plot is taken in a 
characteristic (x, z) plane. The extensions are H in the z direction 
and 8H in the x direction. Due to regions of very small magnetic 
field, the range of variation of magnetic energy density is much 
larger than that of the density being as much as ~ 10 7 . These 
plots are typical of what is found for the turbulence once it has 
reached a statistically steady state. Note that trailing elongated 
features are also apparent in the distribution of the magnetic 
energy density. 



Tn e — Tm 



2YX{P/p) 



(26) 



The angular momentum flow across a line of constant x is 
given by 



T = EWaP/p. 



(27) 



Here W = R 2 for the shearing box simulations or W — r 2 
for the global simulations. 



We now present numerical results for the evolution of the 
simulations. 



5.1 Shearing Box Simulations 

All of the shearing boxes are derived from the simulation 
BaO which was carried out without a protoplanet. It settled 
into a quasi steady state with conserved zero net flux and 
dynamo maintained turbulence as has been been described 
by others (e.g. Hawley, Gammie & Balbus 1996). This model 
was run for 100 orbits at the box centre. The ratio of the 
total magnetic energy to the volume integrated pressure, 
|B 2 /(87r)dV/ J PdV = is plotted as a function of 

time in figure 1. 

After running for 354 time units, this model was used 
to provide initial conditions for models with protoplanets of 
varying masses. Simulations Bal, Ba2, Ba3, and Ba4 were 
begun in this way with M P R 3 /(M*H 3 ) = 0.1,0.3,1.0 and 
2.0 respectively. We also plot J B 2 /(87r)dV/ / PdV = (1//3) 
as a function of time in figure 1 for models Bal, Ba2, Ba3, 
and Ba4. The plots begin to separate from each other shortly 
after after the simulations with protoplanets have begun. 
Note that although (1//3) is in general confined between 
0.008 and 0.02 with or without a protoplanet, there are 
strong fluctuations on a range of timescales down to the 
orbital timescale. The existence of such fluctuations and the 
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Figure 6. Density contour plots in a typical (x,y) plane for 
simulations Bal, top left, Ba2 top right, Ba3 bottom left, and 
Ba4 bottom right taken near the end of the simulations. These 
had M P R 3 /(M*H 3 ) = 0.1,0.3,1.0 and 2.0 respectively. As 
M P R 3 1 (M*H 3 ) increases, the wake becomes more prominent, 
material is accreted by the protoplanet and is pushed toward the 
radial boundaries as a gap is formed. The location of the charac- 
teristic ray given by equation (19) is also plotted. The qualitative 
structure seen in these simulations is maintained once a quasi 
steady state has been attained. 



difficulties they cause for denning mean state variables for 
the turbulent flow has been emphasised recently by several 
authors (e.g. papers I and II; Winters, Balbus & Hawley 
2003b; Armitage 2002). 

The volume averaged stress parameter (a) is plotted as 
a function of time for model BaO in figure 2. As for (1//3), 
erratic fluctuations are seen. The sum of the magnetic and 
Reynolds' stresses nearly always exceeds the magnetic con- 
tribution for which the short term fluctuations are signifi- 
cantly less severe. In these respects the behaviour of local 
and global simulations are similar (Steinacker & Papaloizou 
2002, papers I and II, and section 5.2). The running time av- 
erage of (a) is also plotted and this shows much smoother be- 
haviour. Somewhat large initial values occur because of the 



Figure 7. As in figure 6 but in this case magnetic energy den- 
sity contours are plotted. As M p i? 3 /(M*/f 3 ) increases and a gap 
forms, the magnetic energy generally decreases in the gap region. 
However, some remains concentrated in regions near the high den- 
sity wakes. The qualitative structure seen in these simulations is 
maintained once a quasi steady state has been attained. 

initial strong instability which produced a channel phase. 
However, after about four orbits at the centre of the box, 
the running time average is within a factor of two of its final 
value of 0.008. This value is in agreement with that given in 
Winters, Balbus & Hawley (2003b). But note that the run- 
ning time average still shows variations at the 10 percent 
level after 100 orbits. 

The pattern of turbulence in simulation BaO is estab- 
lished after several orbits. A typical density contour plot in 
a (a;, y) plane is presented in figure 3 near the end of the sim- 
ulation. The range of variation in the density in this plane 
is about a factor of three. The situation illustrated here is 
typical of what is obtained in any (x, y) plane once a statis- 
tically steady state turbulence is set up. Elongated trailing 
structures are noticeable. We plot contours of the magnetic 
energy density in figure 4 which correspond to the plot in 
figure 3. The upper plot corresponds to the (x,y) plane il- 
lustrated in figure 3 while the lower plot is taken in a char- 
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strain*. 




Figure 8. Plots corresponding to those in figure 6 but for models 
Bbl, Bb2, Bb3, Bb4. A similar sequence is indicated. 



Figure 9. Plots corresponding to those in figure 7 but for models 
Bbl, Bb2, Bb3, Bb4. A similar sequence is indicated. 



acteristic (x, z) plane. Due to the existence of regions where 
the magnetic field is near zero, the range of the magnetic 
energy density is much larger than that of the mass density, 
being as much as ~ 10 7 . The plots are typical of what is 
found once the turbulence has reached a quasi steady state. 

When a protoplanet is placed in the centre of the box 
waves are excited that cause outward angular momentum 
transport. These waves also have an associated Reynolds' 
stress (e.g. Papaloizou & Lin 1984) and may affect the un- 
derlying turbulence. 

The volume averaged total stress parameter (a) is plot- 
ted as a function of time for models Bal, Ba2, Ba3, and Ba4 
in figure 5. The plots commence from the moment when 
the protoplanet was introduced. The introduction causes 
the excitation of a wave and an increase in the Reynolds' 
stress which in turn increases the total stress. This is par- 
ticularly noticeable in simulations Ba3 and Ba4 which have 
M P R 3 /(M*H 3 ) = 1.0 and 2.0 respectively and so have large 
enough protoplanets to form a gap. This behaviour was also 
found in the simulation of a 5 Jupiter mass planet in a disc 



with H/R = 0.1 presented in paper II and in the global sim- 
ulations described in section 5.2. After a sharp initial rise 
the volume averaged Reynolds' stress here tends to decrease 
as a gap is opened and material is moved away from the 
protoplanet. The much less noisy volume averaged magnetic 
stresses are also plotted in figure 5. These behave similarly 
for each of the simulations Bal, Ba2, Ba3 and Ba4 but they 
then separate at about time 500. Although simulation Ba4 
which has the highest protoplanet mass also tends to have 
the highest volume averaged magnetic stress, there is no 
other such correlation between magnetic stress and proto- 
planet mass and so this may not be significant. 

Although, as we have indicated above, strong fluctu- 
ations are indeed the rule, there are definite behavioural 
trends as the protoplanet mass increases. To illustrate these, 
density contour plots in a typical (x, y) plane for simulations 
Bal, Ba2 Ba3 and Ba4 near the end of the simulations are 
shown in figure 6. It can be seen that as M P R 3 /(M*H 3 ) in- 
creases, the wake becomes more prominent. It can be seen 
that this wake approximately coincides with that given by 
equation (19) derived from the discussion about character- 
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Vertically and azimuthally averaged densit 



0.0015 - 



0.0005 - 




Vertically and azimuthally averaged stress parameter 



Figure 10. Density averaged over y and z plotted against x 
for simulations BaO full curve, Bal dotted curve, Ba2 dashed 
curve, Ba3 dot dashed curve, and Ba4 triple dot dashed 
curve, taken near the end of the simulations. These had 
M P R 3 /(M,H 3 ) = 0.1,0.3,1.0 and 2.0 respectively. The plots 
show that M p R 3 /(M^H 3 ) increases, material is accreted by the 
protoplanct and is pushed toward the radial boundaries as a gap 
is formed. Such a gap is noticeable in simulations Ba3 and Ba4 
but not in simulations Bal and Ba2 which have a very similar 
behaviour of the mean density to that found in simulation BaO 
which had no protoplanet. Although there are turbulent fluctu- 
ations the same structure persists once a quasi steady state is 
achieved. 



istic rays. The wake is the first observable manifestation of 
the presence of the protoplanet and can be just traced in 
simulation Bal for which M p R 3 /(M f H 3 ) = 0.1. This is well 
before any effects can be seen in the magnetic energy density 
contours or a gap forms. In figure 7 we plot the magnetic 
energy density contours corresponding to the density plots 
in figure 6. 

As M P R 3 /(M,H 3 ) increases and a gap forms, the mag- 
netic energy generally decreases in the gap region. Where 
it remains it tends to be concentrated in regions near the 
high density wakes. This tendency was also apparent in the 




Figure 11. The stress parameter averaged over y and z plotted 
against x for simulations BaO full curve, Bal dotted curve, Ba2 
dashed curve, Ba3 dot dashed curve, and Ba4 triple dot dashed 
curve, taken near the end of the simulations. Although there are 
erratic fluctuations, the variation in (a) among the different sim- 
ulations is typically less than that found in the density in the gap 
region. This is consistent with a general reduction of magnetic 
energy density and stress in that region when a gap is formed. 



global simulation of paper II and is also observed in the 
global simulations G4 and G5 described in section 5.2. 

To examine the effect of lengthening the box in the az- 
imuthal or y direction, we carried out simulations Bbl, Bb2, 
Bb3 and Bb4 which had the same resolution as BaO, Bal, 
Ba2, Ba3 and Ba4 but were extended by a factor of 2 to 87r in 
y. However, because of the obvious increased computational 
demands, these simulations were run for shorter times. Each 
was begun by taking simulation BaO at a time 216, extend- 
ing it to 8n in y by using the appropriate periodicity condi- 
tion and inserting protoplanets with masses given in table 1. 
Mass density and magnetic energy density plots are plotted 
in figures 8 and 9 respectively. A similar behavioural trend 
to that found for the shorter boxes is noted. 

To illustrate the gap formation process the density av- 
eraged over y and z is plotted against x for simulations BaO, 
Bal, Ba2, Ba3 and Ba4 at a moment near the end of the 
simulations in figure 10. As M P R j (M,H ) increases, mate- 
rial is pushed toward the radial boundaries by the action of 
the wave stresses and a gap is formed. In addition it is ap- 
parent that there is some accretion onto the protoplanet in 
the sense that some material settles onto it becoming bound. 
The horizontally and vertically averaged density profiles as- 
sociated with BaO, Bal and Ba2 indicate no sign of a gap 
with the variations between them similar to what would be 
seen for any one of them at different times. However a per- 
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Figure 12. The density averaged over y and z plotted against 
x for simulation Ba4 full curve, the laminar disc simulation with 
a = 0.008, dashed curve and the inviscid laminar simulation, 
dot dashed curve, taken near the end of the simulations. The 
inviscid case has the deepest gap and the case with a = 0.008 the 
narrowest. 



sistent gap is clearly visible for runs Ba3 and Ba4 for which 
M P R 3 /(M*H 3 ) > 1. 

The horizontally and vertically averaged stress parame- 
ter corresponding to the plots in figure 10 is plotted against 
x for simulations BaO, Bal, Ba2, Ba3 and Ba4 in figure 11. 
Although there are erratic fluctuations, the range of vari- 
ation in (a) between the simulations is typically less than 
that found in the density when there is a gap. This is con- 
sistent with a tendency for there to be a general reduction 
of magnetic energy and stress in a gap region that roughly 
scales with the decrease in gas density and pressure there. 

We have also compared the structure of the gap near the 
end of the simulation Ba4 with that obtained from a laminar 
simulation with no magnetic field but uniform anomalous 
Navier Stokes viscosity corresponding to the time average 
value of a = 0.008 obtained from run BaO which had no 
planet. We emphasize that such a time average has no clear 
physical connection with the behaviour of simulation Ba4 
nonetheless we perform the comparison for general interest. 

In figure 12 we plot the azimuthally and vertically av- 
eraged density profile at the end of simulation Ba4 together 
with plots obtained from laminar disc runs with a — 0.008 
and no viscosity which were run for the same times with the 
perturbing planet which had M P R 3 /(M*H 3 ) = 2. These 
have been scaled to have the same magnitude at the box 
edges. 

As expected, the inviscid non magnetic run produces 
a deeper and wider gap than the others. The run with the 
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Figure 13. Velocity vectors in a typical (x, y) plane near the end 
of simulation Bal with M P R 3 /(M*H 3 ) = 0.1. The perturbation 
due to the protoplanet is too weak for the existence of horse shoe 
trajectories that do not cross the whole y domain to be manifest. 

imposed a = 0.008 produces a somewhat narrower gap than 
the MHD run Ba4 which in this way behaves as if it had 
lower viscosity. This effect was also found for the global sim- 
ulation with 5 Jupiter mases in paper II. However, in our 
case there is slightly more material accreted onto the planet 
than the laminar viscous case. 

Finally we examine the flow in the coorbital region 
to check for the appearance of circulating horse shoe 
like trajectories in that region. Velocity vectors in typical 
(x, y) planes for simulations Bal, weakly perturbed with 
M P R 3 /(M*H 3 } = 0.1 and Ba3, strongly perturbed with 
M P R 3 /(M*H 3 ) = 2 are presented in figures 13 and 14 re- 
spectively. Circulating trajectories are apparent in the case 
of Ba3 indicating the existence of a coorbital dynamics. 
However, the flow pattern appears complex with circulating 
eddies and multiple shocks. Furthermore we find no indi- 
cation of a circulating flow around the protoplanet within 
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Figure 14. Velocity vectors in a typical (x, y) plane for simula- 
tion Ba3 with M P R 3 /(M,H 3 ) = 2. The perturbation due to the 
protoplanct is sufficient for the existence of horse shoe trajectories 
restricted to a subset of the y domain to be manifest. 

the Hill sphere. This may be because of a magnetic braking 
effect already discussed in paper II (also see discussion in 
section 5.2). 

5.2 Global Simulations 

In this section we describe the results obtained from the 
global simulations. We begin by discussing the underlying 
turbulent disc model, before describing the interaction be- 
tween protoplanets and global turbulent discs. 

5.2.1 Global Turbulent Disc Model 

The physical parameters of the global disc model are de- 
scribed in section 2.1. The initial disc model had a 7r/4 az- 
imuthal domain, with an imposed zero-net flux toroidal field 
that varied sinusoidally with radius. The initial value of the 



Figure 16. This figure shows the evolution of (a) for the global 
disc model. The solid line represents the Maxwell stress, the dot- 
ted line represents the Reynolds stress, and the dashed line rep- 
resents the total stress. 

ratio of the total magnetic energy to the volume integrated 
pressure, /B 2 /(87r)dV/ / PdV = (1//3) = 0.032 The time 
evolution of this quantity is shown in figure 15. It is clear 
that the initially high magnetic energy decreases as the ini- 
tially imposed magnetic field undergoes reconnection, and 
(1//3) reaches a saturated value of (1//3) ^ 0.011 once the 
turbulence reaches a statistical steady state. 

Figure 16 shows the time evolution of the volume aver- 
aged Shakura-Sunyaev stress parameter a defined by equa- 
tion 26. The solid line represents the Maxwell stress, the 
dotted line represents the Reynolds stress, and the dashed 
line represents the sum of these. It is clear that as the disc 
becomes unstable to the MRI, the stresses increase initially, 
but then decrease to a saturated state corresponding to 
a ~ 7 x 10" 3 . 

The radial distribution of a, time averaged for a short 
period of time (one orbit at r = 1) is shown in figure 17. As 
expected this is a rapidly varying quantity as a function of 
both space and time due to the turbulence (see also paper 
I)- 

The final state of the disc model shown in figures 15-17 
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Figure 17. This figure shows the radial profile of (a) for the 
global disc model. The heavy solid line represents the Maxwell 
stress, the dotted line represents the Reynolds stress, and the 
thin solid line represents the total stress. 

was used as the initial condition for simulations that exam- 
ine the interaction between turbulent discs and embedded 
protoplanets. Models Gl - G4 used full 2-k azimuthal do- 
mains that were constructed by patching eight copies of the 
7r/4 turbulent disc model together. Model G5, which had a 
restricted azimuthal domain of tt/2, consisted of two copies 
of the relaxed, turbulent disc model patched together, with 
periodic boundary conditions in the azimuthal direction be- 
ing employed. These simulations are described below. 

5.2.2 Global Models with Embedded Protoplanets 

In this section we present the results from simulations la- 
beled as Gl - G5 in table 2. The aim is to examine how the 
disc properties change as a function of planet mass (or alter- 
natively as a function of (M p /M t ) /(H/r) 3 ), moving from a 
situation where the protoplanet induces linear perturbations 
in the disc (e.g. run Gl) through to a scenario in which the 
planet perturbation is nonlinear (e.g. run G5). In particular, 
we are interested in how the appearance of the spiral den- 
sity waves vary, how the morphology of the magnetic field 
changes in the vicinity of the planet, when gap formation 
arises in a global disc model, and the impact of the proto- 
planet on the velocity field in the disc. We are also interested 
in how the results of similar local and global simulations 
compare. 

When discussing the results of the global simulations 
we will concentrate on the runs Gl, G2, G3, and G5. Run 
G4 adopted the exact same physical parameters as run G5, 
except that the azimuthal domain in G5 ran from to 7r/2, 
rather than being a full 2n domain. This means that we were 
able to run simulation G5 for ~ 60 planetary orbits because 
of the lower computational cost, whereas run G4 was only 
run for ~ 11 planetary orbits. We therefore give the results 
of run G5 higher prominence in our discussion, whilst noting 
that the results obtained in run G4 are in good qualitative 
agreement with G5. 

Figure 18 shows contour plots of the midplane density 
distribution in the vicinity of the protoplanet for the runs 
Gl, G2, G3, and G5, respectively. In the first panel it is clear 



that the presence of the planet in the disc is essentially unde- 
tectable since the density perturbations induced by the tur- 
bulence are of higher amplitude than the spiral waves excited 
by the protoplanet. The value of (M p /M*)/(H/r) 3 in this 
case is 0.03, with the protoplanet mass corresponding to ~ 3 
Earth masses (assuming a solar mass central star). The sec- 
ond panel, corresponding to run G2, shows that the presence 
of the protoplanet is just about detectable, but the ampli- 
tude of perturbations due to the turbulence still exceed the 
spiral waves induced by the protoplanet. The protoplanet 
mass in this case is equivalent to ~ 10 Earth masses. The 
third panel corresponds to run G3, for which the protoplanet 
mass is equivalent to ~ 30 Earth masses. The presence of the 
protoplanet in this case is clearly detectable, with the spiral 
waves being of similar amplitude to the turbulent wakes. The 
fourth panel shows the results from run G5, in which the pro- 
toplanet mass corresponds to 3 Jupiter masses. In this case 
the disc-protoplanet interaction leads to gap formation since 
the interaction is strongly nonlinear ((M p /M*)/(H/r) 3 ) ~ 8 
in this case) . Overall the results here are in very good agree- 
ment with those obtained in the shearing box runs, in that 
they show the same consistent pattern of behaviour with 
increasing values of (M p /Mt,)/(H/r) 3 . 

We next consider the effect of the protoplanet on the 
magnetic field. Figure 19 shows contour plots of |B| 2 for the 
runs Gl, G2, G3, and G5. Strong variation in magnetic field 
strength arises due to the turbulence, as found in the shear- 
ing box runs, and the plots show that the magnetic field is 
largely unaffected by the presence of the embedded proto- 
planets in runs Gl, G2, G3. However, the last panel, corre- 
sponding to run G5, clearly shows that the magnetic field 
strength is increased near the spiral shock waves induced by 
the more massive planet, in agreement with results obtained 
in paper II and those obtained by the shearing box runs. In 
general, the magnetic energy and stress decrease in the gap 
region as the density and pressure decrease, but the disor- 
dered, turbulent magnetic field close to the planet becomes 
ordered and compressed in the spiral wakes, leading to a 
significant increase in the magnetic stress there. This effect 
is shown more clearly in figure 20 which shows a snapshot 
of the magnetic field vectors in run G5. This illustrates how 
the field becomes compressed and ordered as the gas travels 
through the spiral shocks, and also shows how the magnetic 
field is advected into the protoplanet Hill sphere with gas 
that accretes onto the protoplanet. The effect that the field 
modification has on the vertically and azimuthally averaged 
magnetic stress is illustrated by figure 21 which shows the 
time averaged contribution to a from the magnetic stress as 
a function of radius. It is clear that the compression of the 
field along with reduction of the density and pressure in the 
gap leads to an enhanced a value there. 

We now consider how the presence of the protoplanet 
affects the radial density distribution in the disc. Figure 22 
shows the azimuthally averaged midplane density distribu- 
tion for the runs Gl, G2, G3, and G5, as well as the ini- 
tial 1/r density profile. This plot shows that some accretion 
though the disc has occurred during its initial relaxation 
to a final turbulent state. The lower mass protoplanets in 
runs Gl, G2, and G3 have essentially no effect on the un- 
derlying density structure, but a deep gap is formed in run 
G5. This is in line with the expectations discussed in sec- 
tion 2.3 since (M p / 'M,) / \H /rf < 1 for Gl, G2, and G3, 
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Figure 18. This figure shows the midplane density distribution of runs Gl, G2, G3, G5, respectively. The progression of increasing 
planet mass leads to a clear trend in the disturbance experienced by the disc, with the spiral waves becoming increasingly apparent. Run 
G5 leads to the formation of a clear gap. 
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Figure 19. This figure shows the distribution of \B\ 2 in the disc midplane for the runs Gl, G2, G3, G5, respectively. For runs Gl, G2, 
G3 the protoplanct make no discernible perturbation to the magnetic field structure. For run G5 he magnetic field strength is increased 
in the wakes generated by the planet. 
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Density versus Radius 



Nlp-3 klJ - Turb - TTms-1 543.147 





Figure 22. This figure shows the radial density distribution for 
the runs Gl, G2, G3, and G5. The solid line shows the initial 
l/r density profile of the unperturbed disc. The dotted line cor- 
responds to run Gl, the dashed line to G2, and the dash-dotted 
line to G3. It is clear that there is no gap forming in these cases. 
The dash-dot-dot-dotted line corresponds to run G5 which shows 
clear gap formation. 



Figure 20. This figure shows magnetic field vectors in the vicin- 
ity of the protoplanet for run G5. The compression and ordering 
of the field due to the spiral shocks is apparent, as is the advection 
of magnetic flux into the planet Hill sphere. 



Time Averaged Magnetic a versus Radius 




Figure 21. This figure shows a time average of the magnetic 
contribution to a for run G5. The time averaging occurred for 18 
planet orbits. The increase in a in the gap region occurs because 
of the field concentration in the spiral wakes described in the text. 



but {Mp/M^KH/r) 3 ~ 8 for run G5, and is in very good 
agreement with the shearing box runs. 

We now turn to the question of how the global ener- 
getics of the turbulence and the stresses arising from it are 
affected by the presence of the protoplanets. In the case of 
the lower mass planets represented by runs Gl, G2, and 



G3 the global magnetic energy of the simulations in units 
of the volume integrated pressure is essentially unaffected. 
This is in agreement with the results obtained for a massive 
5 Jupiter mass planet embedded in a disc with H/r = 0.1 
presented in paper II and the shearing box runs. Figure 23 
illustrates this with a plot of (1//3) versus time for runs Gl, 
G2, and G3. It is clear that the magnetic energy remains 
close to its original starting value, with the level of variation 
no larger than the fluctuations normally associated with the 
turbulent energy. Simulation G3 was initiated with an un- 
derlying disc model that had relaxed for a slightly shorter 
period of time than models Gl or G2, which were initiated 
with identical underlying turbulent disc models. Although 
the qualitative evolution of the global magnetic energy is 
the same in the three cases, it is also interesting to compare 
the detailed evolution of (1//3) for models Gl and G2. Fig- 
ure 23 shows that initially the evolution is almost identical. 
However, after ~ 100 time units they clearly begin to di- 
verge due to the 'chaotic' nature of the turbulence, as noted 
by Winters, Balbus, & Hawley (2003b). A similar result is 
obtained in the shearing box runs. Nonetheless the global 
properties of the disc models remain similar, even though 
the local details change. This also has an impact on the way 
in which the local statistics of the turbulence are affected 
by the protoplanet, which in turn affects the detailed evolu- 
tion of the orbit of the planet. These issues are discussed in 
greater length in the accompanying paper IV. 

In the same way that the magnetic energy is largely 
unaffected by the presence of the low mass protoplanets, 
the turbulent stresses are also unaffected. This is illustrated 
by figure 24 which shows the various contributions to a as 
a function of time for run G2 and should be compared with 
figure 16. Similar plots are obtained for runs Gl and G3. 

The picture is a little different when we consider run 
G5. The time evolution of (1/(3) is shown in figure 25, which 
shows that the evolution of the magnetic energy appears to 
decrease during the simulation, before reaching a minimum 
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Magnetic Energy/P versus Time 
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Figure 23. This figure shows the evolution of the magnetic en- 
ergy in units of the volume integrated pressure for runs Gl, G2, 
and G3. It is clear that the presence of the planets has little effect 
on the turbulent dynamo as the magnetic energy remains close 
to its initial value. Variations occur with amplitudes no greater 
than those associated with the usual turbulent fluctuations. 
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Figure 24. This figure shows the evolution of the stress parame- 
ter a for the run G2. The upper line represents the total stress, the 
middle line the Maxwell stress, and the lowest line the Reynolds 
stress. The total a ~ 7 X 10 -3 . 

value of (1//3) — 0.005 beyond which it decreases no fur- 
ther. Such a decrease did not occur during the other global 
runs (G1-G3), or in the calculation of a massive planet in 
a turbulent disc in paper II. At the present time it is un- 
clear whether the turbulent energy is being affected by the 
protoplanet in such a manner that the turbulent dynamo is 
operating less efficiently due to the ongoing planetary per- 
turbation, or whether the decrease is the result of a large 
but temporary fluctuation induced by inserting the planet 
instantaneously in the turbulent disc. The former scenario is 
conceivable when one considers that the azimuthal domain 
in run G5 is 7r/2 so that the fluid experiences the (strong) 
perturbation of the planet four times more frequently than 
it would do in a full 2n disc. Such a strong and frequent 
perturbation may be able to affect the underlying dynamo 
in such a way as to reduce the magnetic energy. The shear- 
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Figure 25. This plot shows the evolution of the magnetic en- 
ergy, in units of the volume integrated pressure, for run G5. This 
quantity decreases initially from ~ 0.01 to ~ 0.045 before leveling 
off. 



ing box runs in general do not support this picture, but run 
Bb4 did also show a reduction in magnetic energy once the 
planet was inserted. It is interesting to note that the results 
of Winters, Balbus, & Hawley (2003a) also contain a sugges- 
tion that the saturation level of the MHD turbulence may 
be affected by the presence of a giant planet. 

We now consider the impact of the embedded proto- 
planets on the velocity field of the disc, and specifically the 
point at which the horseshoe orbits are clearly established. 
Figure 26 shows the fluid trajectories for two runs. The up- 
per panel corresponds to a global calculation in which a 
protoplanet with Mp/M 4 = 10 -4 is embedded in a laminar 
disc (i.e. the protoplanet is the same as that in run G3). 
Close inspection of this plot shows that the horseshoe tra- 
jectories are established in this case. The lower panel shows 
the velocity field for run G3, where it is apparent that the 
horseshoe trajectories are disrupted by the turbulent veloc- 
ity field. In other words the turbulence determines the fluid 
trajectories in this case, and the gravitational potential due 
to the planet is unable to establish the horseshoe orbits in 
the coorbital zone that are obtained in a laminar disc model. 

Figure 27 shows the fluid trajectories for run G5. The 
horseshoes orbits in this case are very clearly denned, in 
agreement with the shearing box run Ba4. In this model the 
circulating region around the planet in the Hill sphere is 
also clearly visible. The effect of the turbulence on the ve- 
locity field in such a strongly perturbed model is essentially 
indiscernible in the near vicinity of the planet. 

In paper II the circulating region around the proto- 
planet in the Hill sphere was found to be disrupted, and 
this was tentatively ascribed to magnetic braking caused by 
magnetic linkage between the circumplanetary disc and the 
protostellar disc (see also figure 20) . A similar situation was 
also found in the shearing box runs Ba3 and Ba4, where 
the usual circulating region in the Hill sphere was found to 
be absent. In these runs the gravitational softening adopted 
was b ~ 0.3//, which is quite large, and results in gas that 
accretes into the Hill sphere forming an 'atmosphere' that is 
largely pressure supported but with some angular momen- 
tum. In this case the removal of angular momentum will 
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Figure 26. This plot shows the fluid trajectories for a laminar 
run with Mp/M* = 10~ 4 (upper panel) and for run G3 (lower 
panel). The usual horseshoe orbits can be seen in the upper panel, 
but the turbulence appears to dominate the trajectories in the 
lower panel. Note that the protoplanet is located at (r p ,4> p ) = 
(~ 3,7r). We note that plotting the velocity field for run G3 at 
different times leads to significant variation in the appearance of 
the velocity field because of the turbulence, but that the flow 
morphology in the laminar case is essentially independent of time 
once established. 



lead to a reduction of the spin of the atmosphere, as ob- 
served. The softening used in run G5 was b = 0.1H, and so 
the formation of a rotationally supported circumplanetary 
disc is more pronounced. The removal of angular momen- 
tum in this case should allow further gas accretion into the 
Hill sphere without a modification of the rotation profile. An 
MHD simulation should therefore result in more material ac- 
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Figure 27. This plot shows the fluid trajectories in the vicinity 
of the protoplanet for run G5. The horseshoe orbits are clearly 
visible, as is the circulating region within the Hill sphere. 



creting into the planetary Hill sphere than occurs in a non 
magnetised run, if magnetic braking is indeed important. In 
order to test this we performed a simulation that was similar 
to run G5, except that magnetic fields were neglected and 
an a viscosity was employed with a = 7 x 1CP 3 in a laminar 
disc model. This laminar disc calculation was evolved for an 
identical amount of time to run G5 (t = 1543.14 Q" 1 , equiv- 
alent to ~ 62.12 planetary orbits), and resulted in less mass 
being accreted into the planetary Hill sphere than occurred 
for the magnetised run, as shown in figure 28. This suggests 
that magnetic fields in turbulent discs will play a signifcant 
role in determining the gas accretion rate onto forming gas 
giant planets. 



6 DISCUSSION 

In this paper we have performed both global cylindrical disc 
simulations and local shearing box simulations of protoplan- 
ets interacting with a disc undergoing MHD turbulence with 
zero net flux fields. The aim has been to extend the results 
obtained in a previous paper (paper II) to a wider range 
of protoplanet masses and conditions of perturbation of the 
surrounding disc. 

Global simulations are naturally more realistic but their 
computational demands mean that only very few can be car- 
ried out. The advantage of local shearing box calculations is 
that more runs can be done at higher resolution even though 
larger boxes than are normally considered (e.g. Brandenburg 
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Average Density versus Radius 
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Figure 28. This plot shows the azimuthally averaged midplane 
density as a function of radius for run G5 (dashed line), and for a 
similar run with a viscous but laminar disc in which a = 7 X 10~ 3 
(solid line). The plots correspond to a run time of t = 1543.14. It 
is clear that more gas has accumulated within the Hill sphere of 
run G5 than for the laminar disc run, suggesting that magnetic 
braking of the circumplanetary disc is playing a significant role 
in modulating the accretion rate onto the planet. 



come apparent as the gravity of the protoplanet dominates 
the turbulence. A region in which fluid circulates about the 
protoplanet may be formed in some cases and this may con- 
tain a magnetic field that circulates about the protoplanet 
and links to the wakes, leading to the possibility of mag- 
netic braking of the circumplanetary disc. In this situation 
the two sides of the disc tend to separate as the gap forms, 
signaling the onset of type II migration (e.g. Ward 1997). 

The nature of the migration induced by disc interac- 
tions is of key importance during the later stages of planet 
formation. This has yet to be considered for a disc with 
MHD turbulence for low mass embedded protoplanets. The 
simulation presented in paper II indicates that migration at 
the expected rate occurs when a gap is fully formed in a 
turbulent disc, in agreement with type II migration theory. 
We comment here that the migration rate induced by the 
interaction is a complex issue on account of the strong tur- 
bulent fluctuations, especially for embedded protoplanets. 
This issue is examined in a companion paper (paper IV). 
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et al. 1995; Hawley et al. 1995) are required in order to ad- 
equately contain the response to a perturbing protoplanet. 

Another advantage is that for zero net flux fields 
there exists a natural scaling indicating that results de- 
pend only on the parameters M' p = M P R 3 /(M*H 3 ) and 
Y/H. Using simple dimensional considerations of the con- 
dition for non linear response and balance of viscous and 
tidal torques we estimated a simple condition for gap for- 
mation given by equation 13 in the form M P R 3 /(M*H 3 ) > 
max(Ct, 0.07(Y/H)). For the conditions of the local and 
global simulations (Y = irR) we have performed both here 
and in paper II, this always results in a condition that 
M P R 3 / (Mi,H 3 ) exceed a number of order unity (i.e. the ther- 
mal condition mentioned by Lin & Papaloizou 1993). 

The pattern of behaviour we have found from our sim- 
ulations is the same for both local and global simulations 
(which show excellent agreement in the results obtained 
for the local interaction between protoplanet and turbu- 
lent disc) and it agrees with the notion described above. 
As M P R 3 /(M^H 3 ) is increased to ~ 0.1 the presence of the 
protoplanet is first manifest with the appearance of the well 
known trailing wake that can be identified with a character- 
istic ray emanating from the protoplanet. This wake appears 
to be well defined even though its form may vary erratically 
due to the turbulent fluctuations. At this stage the magnetic 
field in the disc is relatively undisturbed by the planet, and 
the usual horseshoe orbits in the coorbital region are poorly 
defined because the fluid trajectories are still dominated by 
the turbulence. 

When M P R 3 /(M„H 3 ) exceeds a number around unity 
a gap starts to develop inside which the magnetic energy 
density is on average less than that in the rest of the fluid. 
However, in the gap region it becomes concentrated in the 
high density wakes that are generated by the protoplanet. 
At this point the horseshoe orbits in the coorbital zone be- 
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